{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Using Quadratic Unconstrained Binary Optimasation (QUBO) to solve number partitioning\n",
    "## Introduction\n",
    "As has been shown in the presentation, the 'energy' of a system can be\n",
    "_officially_ defined as:\n",
    "\n",
    "\\begin{equation}\n",
    "    E = \\frac{1}{2}Q \\cdot (\\boldsymbol{x} \\times \\boldsymbol{x}) +\n",
    "     l \\cdot \\boldsymbol{x}.\n",
    "\\end{equation}\n",
    "\n",
    "However, since we choose to depart from the scientifically accurate Ising model\n",
    "toward non-physical problems, we can reduce this equation and simplify it to\n",
    "prevent confusion. If we use the fact that the $x_i \\in [0, 1]$, then we use\n",
    "that ${x_i}^2 = x_i$, and we can use the _diagonal_ of the matrix $Q$ to save\n",
    "our linear terms $l$. The equation becomes much shorter if we write it as:\n",
    "\n",
    "\\begin{equation}\n",
    "    y = \\sum_{i, j}^n Q_{ij} x_i x_j.\n",
    "\\end{equation}\n",
    "\n",
    "Since the QUBO minimises $y$, it is now the goal to translate a problem into a\n",
    "$Q$ matrix such that the desired result is found in the minimum of $y$. The \n",
    "annealer will handle all the increasing and decreasing in energy for us, \n",
    "so we can focus on only the $Q$.\n",
    "\n",
    "This notebook will handle the case of number partitioning, which will be\n",
    "explained in the following section.\n",
    "\n",
    "## Number partitioning\n",
    "### Description\n",
    "Suppose we have a list of numbers $\\mathcal{S} = \\{s_1, s_2, \\ldots, s_n\\}$, how can we\n",
    "divide the list into two subsets, such that the sum of each subset is\n",
    "equal?\n",
    "\n",
    "### Execution\n",
    "Let $x_i = 0$ if $s_i$ is in group 0, and $1$ if otherwise. Let $A_0$ be the sum\n",
    "of group 0, and $A_1$ the sum of group 1, we can define these as:\n",
    "\n",
    "\\begin{equation}\n",
    "    A_1 = \\sum_i^n x_i s_i, ~ A_0 = \\sum_i^n s_i - x_i s_i.\n",
    "\\end{equation}\n",
    "\n",
    "Since we want both sums to be equal, we want the difference\n",
    "$\\Delta = A_0 - A_1$ to be zero, so this would be a proper value to minimise.\n",
    "However, the minimum of $\\Delta$ is the largest negative number, representing \n",
    "the largest difference. To solve this, we will find the minimum of \n",
    "$\\Delta^2$:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\Delta^2 = \\left( \\sum_i^n s_i - 2 x_i s_i \\right)^2\n",
    "             = \\left( \\mathcal{A} - 2 \\sum_i^n x_i s_i \\right) ^2 ,\n",
    "\\end{equation}\n",
    "\n",
    "where we have written $\\sum_i^n s_i = \\mathcal{A}$ as the full sum of\n",
    "$\\mathcal{S}$. This equation turns out to be writeable as:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\Delta^2 = \\mathcal{A}^2 + 4 Q \\cdot (x \\times x),\n",
    "\\end{equation}\n",
    "\n",
    "where appendix A1 goes into the full derivation to this equation. We only \n",
    "need to know for now that $Q_{ii} = s_i (s_i - \\mathcal{A})$ for the \n",
    "diagonal (or 'linear') terms, and $Q_{ij} = s_i s_j$ for the off-diagonal\n",
    "(or 'quadratic') terms. We can now define $y = Q \\cdot x \\times x$ and \n",
    "use a QUBO-algorithm to find $\\min y$.\n",
    "\n",
    "### Code\n",
    "We start by importing the important packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import dimod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`numpy` will speak for itself. `dimod` is the package which will simulate\n",
    "the annealer for us. We _can_ run our program on a DWave machine, using\n",
    "`dwave-ocean-sdk`, but this requires a [full set-up][1] which we will not go\n",
    "into at the moment.\n",
    "\n",
    "#### Creating an Ising or QUBO model using `dimod`\n",
    "In `dimod`, we can create a so-called `BinaryQuadraticModel` which, according\n",
    "to the `dimod` [docs][2], \"contains Ising and [QUBO] models used by samplers\n",
    "such as the D-Wave system\". We can call each of these models via:\n",
    "\n",
    "[1]: https://docs.ocean.dwavesys.com/en/latest/getting_started.html\n",
    "[2]: https://docs.ocean.dwavesys.com/projects/dimod/en/latest/reference/bqm/binary_quadratic_model.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "ising_model = dimod.BinaryQuadraticModel({}, {}, 0.0, dimod.SPIN)\n",
    "qubo_model = dimod.BinaryQuadraticModel({}, {}, 0.0, dimod.BINARY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The empty `dicts` and the `0.0` passed along in the constructor are the\n",
    "_linear_ terms $Q_{ii}$, _quadratic_ terms $Q_{ij}$ and the _offset_,\n",
    "respectively. We ignore the latter for this case.\n",
    "\n",
    "As we have found in the introduction, our $Q$ looks as follows:\n",
    "\\begin{equation}\n",
    "    Q_{ii} = s_i (s_i - \\mathcal{A}); ~ Q_{ij} = s_i s_j,\n",
    "\\end{equation}\n",
    "so we will translate this into one matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# this method is not optimal, but structured for readability\n",
    "def qubo_partition_matrix(numbers: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"Takes a list of numbers which need to be partitioned and returns the \n",
    "    interaction matrix Q to pass to a QUBO.\"\"\"\n",
    "    matrix = np.zeros((len(numbers), len(numbers)))\n",
    "    linear = numbers * (numbers - np.sum(numbers))  # becomes an n array\n",
    "    quadratic = np.outer(numbers, numbers)  # becomes an nxn array\n",
    "\n",
    "    matrix = quadratic\n",
    "    np.fill_diagonal(matrix, linear)  # overwrite diagonal with linear terms\n",
    "\n",
    "    return matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now define the list of numbers which we want to partition, and create the\n",
    "matrix $Q$, and create our `BinaryQuadraticModel` from this matrix via\n",
    "`from_numpy_matrix()`. This method will automatically return a binary model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ -9   4   5]\n",
      " [  4 -24  20]\n",
      " [  5  20 -25]]\n"
     ]
    }
   ],
   "source": [
    "to_partition = np.array([1, 4, 5])\n",
    "interaction_matrix = qubo_partition_matrix(to_partition)\n",
    "print(interaction_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "model = dimod.BinaryQuadraticModel.from_numpy_matrix(interaction_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "As you can see, we have used a different constructor than in the example \n",
    "above. I personally prefer to define a matrix first, and then pass this \n",
    "on to the constructor, instead of fiddling around with dictionaries.\n",
    "\n",
    "We can now start the sampling. For the sampling, we use the `ExactSolver`,\n",
    "which solves the QUBO classically on your local CPU. If we use this solver\n",
    "together with its `sample()` method, we get a result which is comparable to\n",
    "what D-Wave would give us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   0  1  2 energy num_oc.\n",
      "2  1  1  0  -25.0       1\n",
      "7  0  0  1  -25.0       1\n",
      "3  0  1  0  -24.0       1\n",
      "6  1  0  1  -24.0       1\n",
      "1  1  0  0   -9.0       1\n",
      "4  0  1  1   -9.0       1\n",
      "0  0  0  0    0.0       1\n",
      "5  1  1  1    0.0       1\n",
      "['BINARY', 8 rows, 8 samples, 3 variables]\n"
     ]
    }
   ],
   "source": [
    "sampleset = dimod.ExactSolver().sample(model)\n",
    "print(sampleset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Please try this with different sets of numbers $\\mathcal{S}$, maybe even\n",
    "sets where the partition cannot be done neatly into two equal sums.\n",
    "\n",
    "#### Collecting thoughts\n",
    "What you might think: if this can be run classically, why take the effort of\n",
    "running this on a quantum computer? That though might be correct, however,\n",
    "we _simulate_ the calculation of the QUBO, which takes large amounts of\n",
    "useless computer power which can be done more optimally in a different way on\n",
    "a classical computer. As noted in the `dimod` [docs][1], the `ExactSolver`\n",
    "becomes slow for problems with 18 or more variables.\n",
    "\n",
    "Furthermore, in Quantum Annealing, we have the power of [Quantum Tunneling][2],\n",
    "which allows the material which is being annealed to transfer _through_ a\n",
    "potential hill to a different state, where a classical material would only\n",
    "be able to reach that state after first having increased the temperature to\n",
    "an energy higher than said potential hill to pass over it. See the image below\n",
    "for a clarification.\n",
    "\n",
    "![quantum tunneling](https://qph.fs.quoracdn.net/main-qimg-674c5917220f56bbaa2d611bb8e1c78f.webp)\n",
    "\n",
    "[1]: https://docs.ocean.dwavesys.com/projects/dimod/en/latest/reference/sampler_composites/samplers.html?highlight=ExactSolver#dimod.reference.samplers.exact_solver.ExactSolver\n",
    "[2]: https://en.wikipedia.org/wiki/Quantum_tunnelling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Appendix\n",
    "### A1: derivation of $\\Delta^2$\n",
    "Please remember that $x_i = {x_i}^2$, since $x_i \\in [0, 1]$.\n",
    "\n",
    "\\begin{align}\n",
    "    \\Delta ^2 &= \\left(\n",
    "                    \\mathcal{A} - 2 \\sum_i^n x_i s_i\n",
    "                 \\right) ^2 \\\\\n",
    "              &= \\mathcal{A}^2\n",
    "                - 4 \\mathcal{A} \\sum_i^n x_i s_i\n",
    "                + 4 \\sum_{i, j}^n x_i s_i x_j s_j \\\\\n",
    "              &= \\mathcal{A}^2\n",
    "                + 4 \\left(\n",
    "                    -\\mathcal{A} \\sum_i^n {x_i}^2 s_i\n",
    "                    + \\sum_{i=j}^n x_i s_i x_j s_j\n",
    "                    + \\sum_{i\\neq j}^n x_i x_j s_i s_j\n",
    "                \\right) \\\\\n",
    "              &= \\mathcal{A}^2\n",
    "                + 4 \\left(\n",
    "                    -\\mathcal{A} \\sum_i^n {x_i}^2 s_i\n",
    "                    + \\sum_i^n {x_i}^2 {s_i}^2\n",
    "                    + \\sum_{i\\neq j}^n x_i x_j s_i s_j\n",
    "                \\right) \\\\\n",
    "              &= \\mathcal{A}^2\n",
    "                + 4 \\left(\n",
    "                    \\sum_i^n {x_i}^2 {s_i}^2 - \\mathcal{A} {x_i}^2 s_i\n",
    "                    + \\sum_{i\\neq j}^n x_i x_j s_i s_j\n",
    "                \\right) \\\\\n",
    "              &= \\mathcal{A}^2\n",
    "                + 4 \\left(\n",
    "                    \\sum_i^n {x_i}^2 \\left(\n",
    "                        {s_i}^2 - \\mathcal{A} s_i\n",
    "                    \\right)\n",
    "                    + \\sum_{i\\neq j}^n x_i x_j s_i s_j\n",
    "                \\right) \\\\\n",
    "              &= \\mathcal{A}^2\n",
    "                + 4 \\left(\n",
    "                    \\sum_i^n {x_i}^2 s_i \\left(\n",
    "                        s_i - \\mathcal{A}\n",
    "                    \\right)\n",
    "                    + \\sum_{i\\neq j}^n x_i x_j s_i s_j\n",
    "                \\right) \\\\\n",
    "              &= \\mathcal{A}^2\n",
    "                + 4 \\left(\n",
    "                    \\sum_i^n {x_i}^2 Q_{ii}\n",
    "                    + \\sum_{i\\neq j}^n x_i x_j Q_{ij}\n",
    "                \\right) \\\\\n",
    "              &= \\mathcal{A}^2\n",
    "                + 4 Q \\cdot (x \\times x),\n",
    "\\end{align}\n",
    "\n",
    "where $Q_{ii} = s_i\\left(s_i - \\mathcal{A}\\right)$ and $Q_{ij} = s_i s_j$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
